/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "blockopt.h"

#include "em.h"
#include "fastonebigheader.h"
#include "qfunc.h"

#include <string.h>

float
blockopt (uint64_t                               t,
          float                                  eta,
          float                                  rho,
          uint64_t*                              last_t,
          unsigned int                           n_workers,
          unsigned int                           n_items,
          float*                                 kappa,
          float*                                 logalpha,
          unsigned int                           n_labels,
          float*                                 gamma,
          float*                                 tau,
          float*                                 dtau,
          float*                                 rho_lambda,
          float*                                 loglambda,
          float*                                 dloglambda,
          float*                                 pz,
          const float*                           logpriorz,
          int                                    clamp,
          const Rating*                          ratings,
          unsigned int                           n_ratings,
          const Distribution*                    priorkappa,
          const ParameterizedDistribution*       priorlogalpha,
          float                                  sigma,
          bool                                   test_only)
{
  /* In the typical crowdsourcing setup a worker can only work on 
   * an item once, so don't try to guard against duplicate workers.
   */

  /* Step 1: apply prior sparsely
   * NB: this is just a (cheap!) approximation
   */

  /* NB: Limit rho->1 exists, but I'm lazy. */

  if (! test_only)
    {
      float powt = fastpow (t, 1.0f - rho);

      for (unsigned int i = 0; i < n_ratings; ++i)
        {
          uint64_t s = last_t[ratings[i].worker];
          float pows = fastpow (s, 1.0f - rho);
          float decay = 
            fastexp (eta * (pows - powt) / (n_items * (1.0f - rho) * sigma));

          /* tau and gamma */

          float* tauptr = tau + ratings[i].worker * (n_labels - 1);

          for (unsigned int k = 0; k + 1 < n_labels; ++k)
            {
              float dtautmp = (1.0 - decay) * (gamma[k] - tauptr[k]);

              tauptr[k] += dtautmp;
              gamma[k] += dtautmp / (1.0 + n_workers);
            }

          /* loglambda and rho_lambda */

          float dloglambdatmp = 
            (1.0 - decay) * 
            (*rho_lambda - loglambda[ratings[i].worker]);

          loglambda[ratings[i].worker] += dloglambdatmp;
          *rho_lambda += dloglambdatmp / (1.0 + n_workers);

          last_t[ratings[i].worker] = t;
        }
    }

  /* Step 2: maximize per-item objective wrt logalpha and pz */

  float q = em (*kappa, logalpha, n_labels, tau, loglambda, 
                pz, logpriorz, clamp, ratings, n_ratings, 
                priorlogalpha, 100, 1e-3);

  if (! test_only)
    {
      /* Step 3: gradient update of tau, loglambda, and kappa */

      for (unsigned int i = 0; i < n_ratings; ++i)
        {
          memset (dtau + ratings[i].worker * (n_labels - 1), 
                  0, 
                  (n_labels - 1) * sizeof (*dtau));
          dloglambda[ratings[i].worker] = 0;
        }

      dqfunc (*logalpha, n_labels, tau, dtau,
              loglambda, dloglambda, pz, ratings, n_ratings);

      float etat = eta * fastpow (t, -rho);

      for (unsigned int i = 0; i < n_ratings; ++i)
        {
          float* tauptr = tau + ratings[i].worker * (n_labels - 1);
          float* dtauptr = dtau + ratings[i].worker * (n_labels - 1);

          /* this is where duplicate workers would cause multiple updates */

          for (unsigned int k = 0; k + 1 < n_labels; ++k)
            {
              tauptr[k] += etat * dtauptr[k];
              gamma[k] += (etat / (1.0 + n_workers)) * dtauptr[k];
            }

          loglambda[ratings[i].worker] += etat * dloglambda[ratings[i].worker];
          *rho_lambda += 
            (etat / (1.0 + n_workers)) * dloglambda[ratings[i].worker];

          *kappa += (etat / n_workers) * 
                    distribution_dlogpdf (priorkappa, *kappa);
          *kappa += (etat / n_workers) * 
                    parameterized_distribution_d1logpdf (priorlogalpha,
                                                         *logalpha,
                                                         *kappa);
        }
    }

  return q;
}

void
finalize (uint64_t                               t,
          float                                  eta,
          float                                  rho,
          uint64_t*                              last_t,
          unsigned int                           n_workers,
          unsigned int                           n_items,
          unsigned int                           n_labels,
          float*                                 gamma,
          float*                                 tau,
          float*                                 rho_lambda,
          float*                                 loglambda,
          float                                  sigma)
{
  /* Step 1: apply prior sparsely
   * NB: this is just a (cheap!) approximation
   */

  /* NB: Limit rho->1 exists, but I'm lazy. */

  float powt = fastpow (t, 1.0f - rho);

  for (unsigned int w = 0; w < n_workers; ++w)
    {
      uint64_t s = last_t[w];
      float pows = fastpow (s, 1.0f - rho);
      float decay = 
        fastexp (eta * (pows - powt) / (n_items * (1.0f - rho) * sigma));

      /* tau and gamma */

      float* tauptr = tau + w * (n_labels - 1);

      for (unsigned int k = 0; k + 1 < n_labels; ++k)
        {
          float dtautmp = (1.0 - decay) * (gamma[k] - tauptr[k]);

          tauptr[k] += dtautmp;
          gamma[k] += dtautmp / (1.0 + n_workers);
        }

      /* loglambda and rho_lambda */

      float dloglambdatmp = 
        (1.0 - decay) * 
        (*rho_lambda - loglambda[w]);

      loglambda[w] += dloglambdatmp;
      *rho_lambda += dloglambdatmp / (1.0 + n_workers);

      last_t[w] = t;
    }
}
